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\ We have investigated the obhquity evolution of terrestrial planets in habitable zones (at ~ 
' lAU) in extrasolar planetary systems, due to tidal interactions with their satellite and host 
^ ' star with wide varieties of satellite-to-planet mass ratio (m/Mp) and initial obliquity (70), 
through numerical calculations and analytical arguments. The obliquity, the angle between 
\0 \ planetary spin axis and its orbit normal, of a terrestrial planet is one of the key factors in 
^ \ determining the planetary surface environments. A recent scenario of terrestrial planet accretion 
i-^' implies that giant impacts of Mars-sized or larger bodies determine the planetary spin and form 
■ ■ satellites. Since the giant impacts would be isotropic, tilted spins (sin 70 ~ 1) are more likely 
to be produced than straight ones (sin 70 ~ 0). The ratio m/Mp is dependent on the impact 
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c/3 . parameters and impactors' mass. However, most of previous studies on tidal evolution of the 
^ ' planet-satellite systems have focused on a particular case of the Earth-Moon systems in which 
m/Mp ~ 0.0125 and 70 ~ 10° or the two-body planar problem in which 70 = 0° and stellar 
^ ' torque is neglected. We numerically integrated the evolution of planetary spin and a satellite 
■ orbit with various m/Mp (from 0.0025 to 0.05) and 70 (from 0° to 180°), taking into account 
the stellar torques and precessional motions of the spin and the orbit. We start with the spin 
axis that almost coincides with the satellite orbit normal, assuming that the spin and the 
satellite are formed by one dominant impact. With initially straight spins, the evolution is 
similar to that of the Earth-Moon system. The satellite monotonically recedes from the planet 
until synchronous state between the spin period and the satellite orbital period is realized. 
The obliquity gradually increases initially but it starts decreasing down to zero as approaching 
the synchronous state. However, we have found that the evolution with initially tiled spins is 
completely different. The satellite's orbit migrates outward with almost constant obliquity until 
the orbit reaches the critical radius ~ 10-20 planetary radii, but then the migration is reversed 
to inward one. At the reversal, the obliquity starts oscillation with large amplitude. The 
oscillation gradually ceases and the obliquity is reduced to ~ 0° during the inward migration. 
The satellite eventually falls onto the planetary surface or it is captured at the synchronous 
state at several planetary radii. We found that the character change of precession about total 
angular momentum vector into that about the planetary orbit normal is responsible for the 
oscillation with large amplitude and the reversal of migration. With the results of numerical 
integration and analytical arguments, we divided the m/Mp-70 space into the regions of the 
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qualitatively different evolution. The peculiar tidal evolution with initially tiled spins give 
deep insights into dynamics of extrasolar planet-satellite systems and discussions of surface 
environments of the planets. 

Key Words: Extrasolar planets - Satellites, dynamics - Celestial mechanics - Rotational 
dynamics - Tide, solid body 

1 Introduction 

Recently, extrasolar planets of 5-20M®, which may be rocky/icy planets, have started being 
discovered by development of radial velocity survey (e.g., Butler et al. 2004, Rivera et al. 
2005, Lovis et al. 2006) and gravitational microlensing survey (Beaulieu et al. 2006), although 
the majority of extrasolar planets so far discovered have masses larger than Saturn's mass 
(> IOOM0), which may be gas giant planets. If the core accretion model (e.g., Mizuno 1980; 
Bodenheimer and Pollack 1986) is responsible for formation of the extrasolar gas giants, their 
high occurrence rate (> 5%) implies the ubiquity of extrasolar terrestrial planets (e.g., Ida and 
Lin 2004), because cores of gas giants are formed through planetesimal accretion and failed 
bodies that are not massive enough for gas accretion onto them are no other than terrestrial 
planets and icy planets. Near-future space telescopes, COROT and KEPLER may find Earth- 
size planets, through transit survey, including those within so-called "habitable zones," where 
planets can maintain liquid water on their surfaces (e.g., Borucki et al. 2003; Ruden 1999). 

In addition to existence of liquid water, stable climate on timescales more than 10^ yrs may 
be one of essential components for planets to be habitable, in particular, for a land-based life. 
Planetary global climate is greatly influenced by insolation distribution (Milankovitch 1941; 
Berger 1984, 1989), which is largely related with obliquity 7, the angle between the spin axis 
and the orbit normal of the planet (Ward 1974). For example, if 7 > 54°, the planet receives 
more annual-averaged insolation at the poles than at the equator and vice versa. For high 7, 
seasonal cycles at high latitudes would also become very pronounced (Williams and Kasting 
1997). Abe et al. (2005) indicated that obliquity is a very important factor for atmospheric 
transport of water to low-latitude areas. 

Planetary obliquity evolves mainly by tidal interactions with its satellite and host star. All of 
the terrestrial planets in our solar system do not maintain their primordial spin state (obliquity 
7 and spin frequency VL). Mercury spins with 7 = and VL of precisely 3/2 times as large as 
its orbital mean motion (Colombo 1965; Colombo and Shapiro 1965). This configuration is 
an outcome of the tidal interaction with the Sun (e.g., Goldreich 1966; Goldreich and Peale 
1966). Venus rotates slowly with a 243-day period and 7 ~ 180°. This spin state could be an 
equilibrium between gravitational and thermal atmospheric tidal torque (Gold and Soter 1969) 
or an outcome of friction at a core- mantle boundary (Goldreich and Peale 1970). Dissipation 
effects combined with planetary perturbations could bring the spin axis to 180° from any initial 
7 (Nerson de Surgy 1996; Yorder 1997; Correia and Laskar 2001). 

The Earth's obliquity is gradually increasing with the receding of the Moon as a consequence 
of the tidal dissipation in Earth mainly induced by the Moon (e.g., Darwin 1879; Goldreich 
1966). Figures [1] show the evolutionary path of the Earth's obliquity (7) and orbital inclina- 
tions of the Moon's orbit to the ecliptic {i) and the Earth's equator (e), which is obtained 
by integrating the present Earth-Moon system back into the past (for integration method, see 
section 2). The ranges of oscillation during precession are indicated by shaded regions. This 
plots will be refereed to later. 
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[Figure [T] 



Mars' obliquity would be suffering from a large-scale oscillation of ~ 25° ±10° on a timescale 
~ 10^-10^ years (Ward 1973, 1974, 1979), by a resonance between a spin precession rate and 
one of the eigenfrequencies of its orbital precession. The spin-orbit resonance is a different 
mechanism to alter planetary spin state, from the tidal evolution. Because Earth's spin preces- 
sion is accelerated out of the resonance by the Moon, Earth's current obliquity fluctuates with 
only ±1.3° around 23.3° (Ward 1974; Laskar and Robutel 1993; Laskar 1996). 

Here, we focus on evolution of planet's spin state and its satellite orbit due to the tidal 
dissipation in the planet caused by the satellite and the host star. The spin-orbit resonance in 
extrasolar planetary systems was addressed in detail elsewhere (Atobe et al. 2004) and will be 
commented on in section 4. When the tide raising body (star and/or satellite) orbits around 
the spinning planet, the planet is deformed by the tide with a certain time interval, resulting 
in a lag angle 5 as illustrated in Fig. O (5 can be either positive or negative.) The attraction 
of the tide raising body yields a torque on the planet and an equal but opposing torque on 
the body. In the case of Fig. [21 the torque retards the planet spin and increases the orbital 
angular momentum and hence semi-major axis of the tide raising body. If the equatorial plane 
of the planet does not coincide with the orbital plane of the tide raising body, this angular 
momentum exchange also changes the obliquity of the planet and the orbital inclination of the 
tide raising body. The tidal torque depends on the mass of the tide raising body and strongly 
on the separation between the bodies. The tide raised on Mercury and Venus are due to the 
Sun, while that on Earth is mainly induced by the Moon rather than the Sun. 

[Figure [2] 

As we will show in later sections, obliquity evolution at ~ lAU is regulated by a satellite 
if the satellite-to-planet mass ratio is > 0.01. The Moon would have been formed by a grazing 
collision with a Mars- mass object during the late stage of Earth accretion (e.g., Stevenson 
1987, Canup 2004). Oligarchic growth model (Kokubo and Ida 1998, 2000) predicts formation 
of isolated Mars-mass bodies at ~ lAU in the case of the minimum mass solar nebula (Hayashi 
1981). The isolated bodies would start orbit crossing by long term distant perturbations on 
timescales longer than than Myrs (Chambers et al. 1996, Iwasaki et al. 2002). Thus, giant 
impacts with objects of more than Mars-mass would be common and satellites may be produced 
if the impacts are grazing ones. The satellite mass is determined by total mass and angular 
momentum of the impact-debris disk (Ida et al. 1997; Kokubo et al. 2000), which would be 
regulated by the impact parameter of the collision and the impactor's mass. 

If planetary spin is mostly determined by the giant impact forming a satellite (Lissauer 
et al. 2000), the spin axis would align with the orbit normal of the satellite. Recent N-body 
simulations of the planet accretion show that satellite forming impacts are almost isotropic 
(Agnor et al. 1999; Chambers 2001). Hence it is expected that the primordial obliquity 70 
has the differential distribution, ^(70)^70 = | sin 70^70, so that the most common initial spins 
are tilted ones (70 ~ 90°). As we will show, the obliquity evolution would be quite different in 
cases of tilted initial spins from the familiar evolution in Fig. [H In order to investigate obliquity 
evolution of extrasolar terrestrial planets, evolution with wide ranges of initial obliquity (70 = 
0°-180°) and satellite mass should be studied. Many studies on obliquity evolution have been 
focusing on the Earth-Moon system of 70 ~ 10°. Although Counselman (1973) and Ward and 
Reid (1973) addressed more general characteristics and outcome of tidal evolution, they assumed 
that the mutual inclinations among the satellite's orbit, planet's equator and orbit are always 
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zero, i.e., 7 = (for details of Counselman (1973), see section 3.2). However, as mentioned 
above, the planetary spin of terrestrial planets is more likely to be initially perpendicular to its 
orbit normal (70 ~ 90°). 

In this study, we investigate the tidal evolution of planet-satellite systems with wide ranges 
of satellite-to-planet mass ratio and initial obliquity, numerically calculating the planet's spin 
state and satellite's orbit. We will show that character change of precessional motions during 
tidal evolution plays an important role in producing diversity of the tidal evolution. In Section 
2, we describe basic equations for precession and tidal evolution. The initial conditions are 
also described. Section 3 describes the results of the simulations and general features of tidal 
evolution. Section 4 and 5 are devoted to discussion and conclusion. 



2 Model and the basic equations 
2.1 Model 

We consider a three-body system composed of a host star with mass M^,, a planet with mass 
Mp and physical radius Rp, and a satellite with mass m. For simplicity, we assume that the 
planet-satellite system is in the circular orbit around the host star with the mean motion Hp 
and that the satellite is in the circular orbit around the planet with the mean motion n. 

We adopted the planet-centric frame {X,Y,Z). In this frame, the satellite and the central 
star rotate around the planet with n and Up, respectively. Figure [3] shows the geometry of the 
orbital and equatorial planes of the planet, and the orbital plane of the satellite, s, k, and 
n represent unit vectors in the directions of spin axis of the planet, its orbit normal, and the 
satellite orbit normal, respectively. In this frame, k is a fixed unit vector in the Z-direction. 
(Although the frame with n fixed might be better to discuss the issue of diversity of tidal 
evolution, we adopt the frame that is more intuitive.) We denote the obliquity of the planet 
(the angle between k and s) by 7, the inclination of the satellite orbit to the planet orbit (the 
angle between k and n) by z, and that to the planetary equator (the angle between n and s) 
by e. We assume the planet as an axisymmetric fluid rotating with geometrical axis always 
parallel to its spin angular momentum vector. The planetary spin angular velocity is denoted 
by Q. The satellite spin angular momentum is neglected (Appendix B). 

[Figure E] 



2.2 Basic equations 

The scalar angular momentum of the planetary spin and that of the satellite orbital motion 
are denoted by H and h. Denoting the total precessional torques acting on the planet and 
the satellite by Lp and Lg and the tidal torques by Tp and Ts, the equations of motion of the 
planetary spin and the satellite orbit are 

^ = L^H-T.. (2) 

The precessional torques just rotate the direction of s and n, keeping H and h constant 
(Eqs. (139|) and (l40l) in Appendix A). The tidal torques alter H and h. Since angular momentum 
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is also exchanged with the planetary orbit, {Hs + hn) does not conserve. But, since the planet 
orbital angular momentum is much larger than H and h, we assume that k is invariant. 

We adopt the method by Goldreich (1966) to calculate Eqs. ([1]) and ([2]) because it follows 
precessional motions, although some formalisms (e.g., Mignard 1978, 1979, 1980) are averaged 
over the precession. The precessional motions play a key role in the tidal evolution of the 
system in the case of high initial 7. (For the Earth-Moon system with relatively low initial 
7, predicted tidal evolution is hardly affected by whether precessional motions are taken into 
account or not.) Goldreich's method is a multiple averaged "secular theory," utilizing three 
distinct timescales: orbital periods, precessional ones, and tidal friction timescale. Because tidal 
evolution timescale is usually much longer than precession periods, Tp and Tg are neglected 
in the equations to describe precessional motions. The equations are (Goldreich 1966; also see 
Appendix A) 

ds 

H— ~ Lp = L(s-n)(sxn)+iri(s-k)(sxk), (3) 

at 

h— ~ Ls = -L(s ■ n)(s X n) + ir2(n • k)(n X k), (4) 

at 

where all the quantities are averaged over the orbital periods (the shortest timescales) of the 
star and the satellite around the planet. H and h are constant with time on the precessional 
periods (Appendix A). The first terms in r. h. s. of Eqs. and express the precessional 
torques between the planet and satellite that cause precession around Hs + hn, and the second 
ones the stellar torques that cause the precession around k. The constants L,Ki, and K2 and 
detailed expressions of Eqs. ([3]) and (jlj) for numerical integration are given in Appendix A. In 
addition to H and h, Az = Hx + hy (the Z-component of total angular momentum in the 
planet-satellite system) and x = Kix^ + K2y^ + Lz^ (a kind of total potential energy) are also 
conserved (Appendix A), where x = s ■ k = cos 7, ?/ = n ■ k = cosi, and z = s ■ n = cose. 
The tidal evolution is 

dHs „ dhn „ , , 

where s and n are averaged over precession periods (while those in Eqs. ([3]) and (jl]) are instan- 
taneous ones), and Lp and Lg vanish by the precession averaging. Tp and Tg are numerically 
averaged over the precession periods. We only consider planetary tides induced by the star and 
the satellite (Appendix B). When the tidal torques are included, H, h, Az, Xi ^i-nd a that are 
treated as constant in Eqs. ([3]) and (jl]) change slowly with time (Eqs. HH [511 [521 [54] in Appendix 
A). We adopt the constant time lag model by Mignard (1981) and Touma and Wisdom (1994) 
for Tp and Tg (Appendix C). 



2.3 Initial conditions 

If the initial planetary spin state is determined by the satellite forming impact, resulting obliq- 
uity would have the distribution stated in Introduction section and the initial spin period would 
be a few hours if perfect accretion is assumed (Agnor et al. 1999; Chambers 2001). N-body 
simulations of the accumulation of a satellite in a circumplanetary disk predict that the satellite 
is formed in an orbital plane close to the planet's equatorial plane (e ~ 0) at an orbital radius 
2.6-4.6i?p (Ida et al. 1997; Kokubo et al. 2000). The mass of the formed satellite depends on 
an impact parameter and the mass of a projectile and a target. 

We simulated the tidal evolution and the evolutionary path of the obliquity and the satellite 
orbit with various m/Mp and initial obliquity 70: m/Mp = 0.0025 x j(j = 1,2, ■ ■ -26), 70 = 
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10 X /(/ = 1, 2, ■ ■ ■ 8, 10, ■ ■ ■ 17). In the numerical calculations, the initial semi-major axis of the 
satellite, oq, and the initial rotation period, Dq, were chosen to be 3.8i?p and 5 hours. Although 
the simulation is limited to these oq and Dq, we will present the dependence of the results on 
Oo and Dq through analytical arguments. As shown later, key quantities to regulate diversity of 
tidal evolution are the radius of precession-type change (ocrit) and the outer co-rotation radius 

('^c,out)- Although flc,out 

oc Dq^ (Eq. [27D, cicrit is dependent on Oq and Dq only weakly (cq does 
not affect a^^out either). As long as Dq and differ from the above values within a factor of a 
few, the overall features of the results presented here do not change. Initial inclination of the 
satellite orbit to the planet's equator, eo, is assumed to 1°. We also carried out the calculation 
with eo = 5° and 10° (results are not shown in this paper). Increase in eo results in an only 
slight increase in the asymptotic value of i at large a and corresponding oscillation of e due to 
precession (compare Fig. [T]with Fig. [5]). However, other properties of evolution are similar. If 
evolution in the case of Fig. [T]is continued, 7 and e will start decreasing to zero as approaching 
the CO- rotation radius. On the assumption that satellites are formed by a giant impact, eo 
would not take larger values. Other parameters are assumed that M^, = IMq, = lAU, 
Mp = IM0, a = 0.33, p = 5.5gcm~^, kg = 0.95, k2 = 0.30, and the constant time lag 6t = 11.5 
min. The dependences of the results on these parameters will be also presented. We integrate 
the evolution until the synchronous state is achieved or the satellite orbit decays to 2Rp. In 
the latter case, we regard that the satellite falls onto the planet. 

2.4 Precessional motions 

In order to understand the numerical results of tidal evolution, we briefly summarize charac- 
teristic behaviors of precessional motions. Three different types of precessional motions are 
shown in Figs. HI which were obtained by numerical integration. The trajectories of s and n are 
projected onto the planetary orbital plane (X, Y), in the Earth- Moon system, for (a) a ~ 5R^, 
(b) a ~ 15i?0, and (c) a ~ 60R^ (the present location), respectively. Although these differ- 
ent precessional motions were already discussed in detail by Goldreich (1966) and Touma & 
Wisdom (1994), we present the summary again, because the character changes in precessional 
motion play an important role to produce diversity of obliquity evolution. 

When a ~ 5-R®, the interaction between the Earth and the Moon is much greater than their 
interaction with the Sun, i.e., L Ki, K2 in Eqs. and (jlj). In this case, s (denoted by a solid 
curve and filled squares) and n (denoted by a dashed curve and filled triangles) precess about 
a common axis with the same precession speeds, resulting in constant e, the angle between s 
and n. The common axis precesses slowly about k {{X,Y) = (0,0)) by the stellar torques, as 
indicated in Fig. H^. Hereafter, we call this motion type I precession. When a ~ 15i?0, torques 
on the Moon's orbit due to the Earth and the Sun are nearly equal, on the other hand, those 
on the Earth is mainly due to the Moon, i.e., K2 ^ L > Ki. In this case, n precesses about 
the average position of s and k, while s tends to precess about n and is dragged by the motion 
of n, as illustrated in Fig. |Dd. We call this motion type II precession. In the current state, 
the Sun's torque on the Moon's orbit is much larger than the Earth's torque, while the Moon's 
torque on the Earth is about 2 times larger than the Sun's torque, i.e., L > Ki. In this 

case, n precesses about k maintaining i almost constant. Although s still tends to precesses 
about n, s actually precesses about k because n precesses about k on much shorter period (18 
years) than the precession period (27000 years) of s about n, as illustrated in Fig. We call 
this motion type HI precession. 

[Figure S] 
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We define the critical orbital radius of the satellite, Ocrit, by K2 = L, at which type I 
precession is transformed into type II. From Eqs. (!35l) and (138!) . 



K2 3 GM* f CL \ 3 fnp\'^ ( a 



(6) 



so that ttcrit is given by 



where D is the spin period of the planet (= 2n/fl) at Ccrit- For the Earth, Ocrit — 17-R® 
(Goldreich 1966). Although D can take various values according to tidal evolution from Dq, 
the dependence on D is weak. 

Because the ratio K2/L varies as the fifth power of [a/Rp), the precession of the satellite is 
almost completely dominated by the planetary torque (type I precession) when a < Ocrit and by 
the stellar torque (type III precession) when a > acrit (Goldreich 1966). Type II precession near 
Ocrit results in changes in all of 7, i, and e. In the case of initially tilted spins, this fluctuation is 
so large that outward migration of the satellite orbit is transformed into inward one (see section 
3). 



3 Results of the tidal evolution 
3.1 Diversity of tidal evolution 

Our numerical calculations show qualitatively different three types (A, B, C) of tidal evolution: 

A) the satellite orbit monotonically expands outward until a reaches an outer co-rotation radius, 

B) it expands to ~ Ocrit and turns back onto the planet, and C) the same as B) but it is locked 
at a synchronous state before falling onto the planet. The tidal evolution of the Earth- Moon 
system belongs to type A evolution, so that characteristics of type A evolution have been well 
known. On the other hand, type B and C evolution is newly found by us. The evolution is type 
B or C if an initial spin is tilted, as shown below. Since tilted initial spins are most common, 
as we stated in Introduction section, the newly found type B and C evolution would be more 
common than type A for extrasolar terrestrial planets. 

3.1.1 Tidal evolution with initially straight spins 

Figures [5] show a typical result of type A evolution with m = 0.02Mp and 70 = 10°: (a) 
obliquity 7, (b) i, (c) e, and (d) the semimajor axis of the satellite (a), the estimated inner 
and co-rotation radii (ac,in and ac,out) (see below and section 3.3). 7, i and e are plotted as 
functions of a. The semimajor axis is plotted as a function of time. The evolution in this case 
is qualitatively similar to that of the Earth- Moon system in Figs. [TJ 

In the proximity of the planet, the spin axis of the planet (s) and the orbit normal of the 
satellite (n) precess about the common axis with nearly constant mutual inclination e (type 
I precession). As long as a < Ocnt ~ 15i?p, the satellite orbit migrates outward maintaining 
e ~ (an initial value) and accordingly 7 ~ i. As the satellite approaches Ccrit, n tends to 
precess about the average position of s and k (type II precession). Because of this transition 
of precession, s and n begin to precess at different rates so that e increases. 



7 



When a > Ocrit, s and n precess about k with different speeds (type III precession), i quickly 
decreases to zero, while 7 starts increasing (Fig. [5]). These behaviors are explained as follows. 
From Eqs. (H9|) and (!50|) in Appendix A with tidal torque formula in Appendix C, 

dx _ 3 Gm^Rlh6tiy-xz) 3 GM.^i^^Mt (1 - x^) 

dt ~ 2 H ^ ' 2 al H ^ ^ 

dy _ 3 Gm'Rlk26t jx-yz) ^ 

dt ~ 2 a6 h 
where x = s ■ k = cos 7, y = n ■ k = cost, z = s ■ n = cose, and St is a time lag for distortion 
of the planet (Appendix C). In Eq. ([8]), the first and the second terms express the stellar and 
satellite's tidal torques. In this region, the satellite tide is dominant. Neglecting the second 
term, dx / dt < ii Qz > 2n, because usually y > xz in this regime. Since fl n except for late 
phase near the co-rotation radius, 7 increases toward the asymptotic value of 7, cos~^(2n/r2), 
which is ~ 90° for Qz/n ^ 2. Equation ([9]) includes only a stellar torque. Since usually 
X > yzm. this regime, dy/dt > and the stellar torque decreases i to zero throughout the tidal 
evolution. 

From Eqs. (153|) and (!54l) with tidal torque formula in Appendix C, evolution of the spin 
frequency fl and the satellite semimajor axis a is 

dn 3Gm^Rlk26t 1 



dt 



p 

P -{fi(l+x2)-2npx}, (10) 



"p -'p 



da 3Gm^R\k2bt^a ,^ , , , 

Tt = 2 (11) 

where /p is the planet's moment of inertia, given by aMpR^ (where a < 2/5; for the Earth, a ~ 
0.33). Neglecting the stellar torque (the second term) in Eq. (ITUi) (assuming rn? ja^ ^ M^/Op), 
dfl/dt < if f2 > (2cose/(l + cos^e))n. Thus, as long as > 2n, the spin rate fl decreases 
for any value of e. Eventually ficose becomes smaller than 2n, so that dx/dt > (Eq. [8]) 
and 7 begins to decrease, as shown in Fig. [St. The asymptotic value is cos~^(2n/r2) — > 0° as 
Vtz/n —>■ 2. For Qz/n < 2, 'j approaches 0° regardless of its value. 

This decrease in 7 causes the stable synchronism {Q ~ n) with 7 ~ i ~ e ~ 0. Equation 
(fTTj) indicates the satellite is receding from the planet when Qz/n > 1. The co-rotation radius 
Oc at which Qz/n = 1 depends on the value of z. There are generally two co-rotation radii in 
the prograde case; the synchronous state is stable at outer one a^out while it is unstable at inner 
one flcin (section 3.2). At a > acrit in type A evolution here, z > 0, so that a monotonically 
increases. As a — a^out; the orbital expansion of the satellite terminates {da/dt = 0). As 
discussed in the above, 7 ~ i ~ (and hence, e ~ 0) as approaching ac,out, so that a^out is 
given hjQ = nm this case. 

Note that the weak stellar tidal torque secularly decreases the planetary spin and hence 
flc.out (Eq. [in])- As a result, the satellite orbit slowly decays, being locked at ac_out (Ward and 
Reid 1973), although the decay timescale is longer than 10^° years for a planet at lAU (e.g., 
Goldreich 1966; also see Eq. [3i 



3.1.2 The cases of initially tilted spins 

[Figure [H] 
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Figures [6] show the evolution in the case of m = O.OlMp and 70 = 80°. During type I 
precession up to ~ 18i?p, 7 and i are kept almost constant. At a ~ 18i?p, type II precession 
starts. As a consequence of transition from type I to type II precession, e suddenly starts 
oscillation with a large amplitude (Fig. Et). During type I precession, since eo = 1° and e is 
almost conserved, s and n always point to the same direction as illustrated in the left panel 
of Fig. [71 In the case of type III precession, s and n tend to precess about k with different 
frequencies. If type III precession starts with 70 and ^ ~ ~ when s and n are at the 
same precession phases, while e ~ 270 when s and n are at the opposite precession phases (the 
right panel of Fig. [7]). When precession is transformed to type II precession at a ~ Ocrit, s 
and n tend to precess differently, although not completely. Hence, when a becomes ~ Ocrit, 
the precession tends to have e oscillate in the range up to ~ 270, which is consistent with the 
numerical result in Fig. [6I3. 

[Figure [7] 

Equation (fTTjl shows that the orbital migration is inward when Qz/n = Qcose/n < 1. 
Because of the large fluctuation of e caused by large 70, Qcose/n < 1 during large part of 
precession. The satellite orbit migrates back and forth during a precession period. In the 
case of Fig. [6], the net migration is inward during type II precession. Approaching the planet, 
the interaction between the planet and satellite dominates the precession again, and s and n 
begin to precess about each other (back to type I precession). The amplitude of precessional 
oscillation of e deceases as the satellite orbit comes back to the type I precession regime. The 
asymptotic value of e is larger than 90° in the case of Fig. El satellite's orbit becomes retrograde 
around the planet (in which Oc vanishes). Hence, satellite eventually falls onto the planet. 

In section I3.2.2[ the reversal of migration at a ~ Ocrit is explained in terms of angular 
momentum and energy and the condition for the reversal is presented. In previous studies, 
only the stellar tidal torque has been considered for a mechanism for eventual decay of initially 
prograde satellite orbits. Here we have newly found that planet-satellite tidal interactions can 
reverse the orbital migration in the case of tilted spins without any help of the stellar tidal 
torques. 

3.1.3 The cases of initially moderately tilted spins 

Figures [8] show the evolution in the case of m = 0.04Mp and 70 = 40°. In this case, Ocrit ~ 14i?p. 
The evolution is similar to that in Fig. [6] until a turns back from Ocrit- In this case, however, 
the satellite's orbital and planetary spin periods are locked at stable synchronism. Figure [Hli 
shows the time evolution of a and the co-rotation radii a^in and Oc^out calculated by Lzo = Lzc 
(Eqs. [231 and [25!) . Since the co- rotation radii are defined by i^cose = n, it largely oscillates 
due to the oscillation of e on the precession period, at a ~ Ocrit- In Fig. [Hli, such oscillation 
occurs from 0.4 x 10^ years to 2.2 x 10^ years. Since the oscillation is so violent, ac,out and 
ac,in during this period is omitted in Fig. [Hli. During this period, ac,out changes so rapidly that 
it passes through the satellite orbit without capturing the orbit at the synchronism (t = 0.5- 
1.0 X 10^ years). After the oscillation ceases, a approaches a^out from the outside of ac,out, and 
the satellite orbit is captured at ac.out- The condition for the capture is analytically derived in 
section I3.2.2[ This evolution is also newly found one in the present paper. 

[Figure [S] 
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3.1.4 High m cases 

If the satellite is massive (m > O.OSMp), the planetary spin and the satellite's orbit become 
synchronous before the occurrence of the precession transition. In this case, 7, i, and e are 
almost conserved as the initial values. We will refer to this type of evolution as type A2 
evolution and the evolution in section 3.1.1 as Al. 

3.1.5 Small m cases 

If the satellite is light enough, the stellar tidal torques dominate over the satellite torques. In 
this case, Q becomes smaller than n before the obliquity becomes zero, then the satellite begins 
to decay toward planet very slowly. The subsequent reduction of the planetary spin leads to a 
synchronous state with planetary mean motion {Q = rip) with 7 ~ i ~ e ~ 0. We will refer to 
this as type A3 evolution. 

3.1.6 Retrograde cases 

If an initial planetary spin is determined by a satellite forming impact, retrograde spins (70 > 
90°) have equal probabihty to prograde spins (70 < 90°). Since we assume that s and n align 
initially (eo ~ 0), io > 90° when 70 > 90°. When the satellite tide is dominant, the second 
terms of r.h.s. of Eqs. ([H]) and (ITOl) are negligible. Compared with the prograde and 
y change sign while z has the same sign. Then evolution of fl and a is the same (Eqs. [10] and 
[TTj) . while the equations for x and y, Eqs. ([8]) and ([9]), change sign. So, the obliquity evolution 
is symmetric about 90°. For 70 > 90°, 7 decreases toward 90° until Qz/n < 2. After that the 
obliquity increases to 180°. In type A3 evolution, however, the secular change in x is dominated 
by the second term in r.h.s. of Eq. ([H]). Then, the obliquity eventually decreases toward 0°. 

3.1.7 Parameter dependence 

Figure [9] shows all the results obtained by numerical calculations. Triangles represent A2 
evolution. Type Al evolution with 7 — 0° (for 70 < 90°) or 7 ^ 180° (for 70 > 90°) are 
plotted by open circles. The results with 7 — 0° for 70 > 90° (type A3), are plotted by filled 
squares. For 70 < 90°, the boundary between Al and A3 is not clear in the numerical results, 
so that A3 for 70 < 90° is also plotted by open circles. Type B and C evolution is expressed by 
crosses and filled circles, respectively. The regions of the qualitatively different tidal evolution 
are clearly divided in the m/Mp-70 plane. In section [3.2.2l we analytically derive the boundaries 
of the regions. 

[Figure [9] 

3.2 Angular momentum and energy 

During the tidal evolution of the system, angular momentum (L) is transfered among the planet, 
the satellite and the host star, while the total mechanical energy (E) decreases through tidal 
dissipation in the planet. In this subsection, we first summarize the arguments in terms of L 
and E for the planar problem developed by Counselman (1972), and generalize it to non-planar 
systems to explain diversity of tidal evolution found by our numerical simulation. 
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3.2.1 Tidal evolution of co-planar planet-satellite system 

If the stellar torques are not included and cq = 0, the tidal evolution is a completely planar 
problem. Because the stellar tides are neglected, type 11 and 111 precessions do not exist. Due 
to e = 0, type 1 precession does not exist, either. In the planar case, total angular momentum 
L and total mechanical energy E of the planet-satellite system are given by 

M m 

L = H + h = aM^RlQ + na^, (12) 

E = -aM^RlQ'-^^. (13) 

These equations are deduced to non-dimensional forms, 

L = n + n-\ (14) 

E = Q^-n^, (15) 

where 

Q = {Q/a)k-^, (16) 

n = {n/ay/^k-\ (17) 

L = L{aMpRlk^a)-\ (18) 

E = E{aMpRlk^a^y\ (19) 

In the above, the frequency a and a parameter k are defined by a = {GM-p/ R^Y^'^ = (AnGp/SY^'^ 
and k = (m/aMp)^/^(l + m/Mp)~^/^^. The condition of spin-orbit synchronism (n = fi) is de- 
duced to 

n = h^. (20) 

At synchronous points the gradient of E is normal to a constant-L contour. 

Contours of L and E and the line VL = 'h? are plotted by solid, dashed, and dotted curves 
in the (f2, n) plane in Fig. [TOl The satellite-planet system evolves along a constant-L line 
(corresponding to an initial value, Lq) in the direction of decreasing E, as indicated by arrows 
on solid curves. Note that smaller \h\ corresponds to larger a. Evolution with decreasing 
(increasing) \n\ represents orbital expansion (decay) of the satellite. 

[Figure [ID] 

The lines of constant-L, constant-L^, and Vt = are mutually tangent at [Vt* ,h*) and 
{—Vt*, —h*), where Q* ^ 0.439 and n* ^ 0.760. If |L| is smaller than the critical value L* = 
L[Cl*,n*) = —L{—Cl*,—n*) ~ 1.755, constant-L line and the synchronism line (Eq. [201) never 
cross each other, which means that the satellite monotonically migrates inward and inevitably 
falls onto the planet. If |L| > L*, there are two synchronous points for each L: the outer points 
(a = ctcout) for \Q\ < Q* are stable and the inner ones (a = ac,in) for \Q\ > Q* are unstable, as 
shown in Fig. [TOl 
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3.2.2 Tidal evolution of non-planar planet-satellite system 

In the non-planar cases (e 7^ 0), the total angular momentum is given by 

L = aMMRls + -^^na'n. (21) 

Since the stellar torque is considered here, L changes through exchange with orbital motion of 
the planet. However, the stellar torque is restricted in the planet's orbital plane {XY plane), 
= L ■ k is conserved, where 

M m 

Lz = aMMRl cos 7 H na^ cos i (22) 

P Mp + m ^ ^ 

Note that Lz — L0COS70, since 7 ~ i initially. 

Figures [TT] show the time evolution of normalized L = |L|, Lz, and E (Eqs. [T8l and [T9l) . in 
the cases of (a) m = 0.02Mp and 70 = 10°, (b) m = O.OlMp and 70 = 80° and (c) m = 0.04Mp 
and 7o = 40°. During the evolution, E monotonically decreases throughout the evolution. 
On the other hand, Lz is conserved in all cases, while Lxy oscillates and secularly decreases 
because of the stellar torque. Therefore L asymptotically approaches Lz (— I/0COS70). Using 
this characteristic, type B evolution is explicitly predicted by initial parameters m/Mp and 70, 
as below. 

[Figure [n] 

In Fig. [T2I we plot the numerically obtained time evolution of normalized Q cos e and n 
(Eqs. [T6l and [T71) . in the cases of (a) m = 0.02Mp and 70 = 10° (type A evolution in section 
3.1.1), (b) m = O.OlMp and 70 = 80° (type B in section 3.1.2) and (c) m = 0.04Mp and 
7o = 40° (type C in section 3.1.3). In these figures, we plot ilcose instead of fl in order to 
fix the line representing spin-orbit synchronism in the non-coplanar case, flcose = n^, on the 
plane. Contours of L are plotted for e = 0. Since e is not zero generally, actual L is different 
from L at instantaneous points of trajectories in the figure except for the initial and final stages 
in which e ~ 0. However, these contours provide good guide for the evolution of the trajectories, 
as shown below. With eo ~ 0, initial positions are in upper right region (n > 0, f2 > 0) or lower 
left region (n < 0, < 0) in Fig. [TOl Since evolution is symmetric between the two regions 
unless stellar tidal torque is dominant, we discuss the cases starting with n > and f2 > 
below. 

In the early phase where type I precession is dominated, e ~ and 7 and i are conserved. As 
a consequence, in this phase, the conservation of Lz implies that of L and a trajectory gradually 
moves along a constant-L line corresponding to Lq. When the satellite approaches Ocrit, stellar 
precession torque becomes important, then e begins to oscillates on precession periods (type II 
precession) and L decreases. Because the horizontal axis is Q cos e in Fig. [12], the oscillation 
is apparent in the figure. When 70 is small, the satellite keeps receding from the planet and 
approaches a^out (type A evolution in Sec. 3.1), as shown in Fig. ITIk . However, in the high 70 
case (Fig. THIh). the oscillation amplitude of e is so large that the trajectory strides over a^out- 
The oscillation timescale is comparable to the precession periods of s. Because it may be much 
shorter than tidal evolution timescale, the system cannot be captured at the synchronism. Once 
the trajectory strides cic,out, orbital migration of the satellite is reversed to inward one. In the 
f^cose-n plane in Figs. [T21 leftward/downward evolution is reversed to rightward/upward one 
(see arrows in Fig. [TUI) . 
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Due to damping of Lxy, L decreases to Lz after the large oscillation. If Lz < r-u I.755, L 
becomes smaller than L*, so that the satellite eventually falls onto the planet (type B evolution 
in Sec. 3.1), as shown in Fig. [T2b . If 70 and hence the oscillation amplitude are not large 
enough (equivalently, Lz — I/0COS70 is not small enough) to enter the regions of L < L^,, 
the satellite meets a^out again during the inward migration (type C evolution in Sec. 3.1), as 
shown in Fig. WDc . Since the oscillation of e has been damped, the system is captured at the 
synchronism. Thus, whether type B or C evolution is predicted by initial Lz'- type B for 
Lz < I/* and C for Lz > L*. 

[Figure [12] 



3.3 Co-rotation radius 

The condition to stride a^out separates type B and C evolution from type A evolution. This 
condition may be roughly obtained by Ocrit ^ Q^JTout' where a^^^^ is the outer co-rotation radius 
evaluated with 7 = z = e = 0. As shown in next subsection, this condition is in good agreement 
with the numerical results. 

The expression of a^out is derived by the conservation of Lz- Lz at qq is 

/ Mm \ 
aM^QoRl + jj^, n{ao)ao' cos 70 = «Mp(]oi?p(l + /o) (23) 

^° - (0:33) (s.Sgcm-s) (omMp) (sh^) (3:81;) ' ^^^^ 

where subscript "0" represents initial values. Using the definition of the co-rotation radius (oc), 
fl{ac) cose = n{ac), the angular momentum at Oc is 



Lzc = C(Mr,n{ac)Rr, \- iTF^, n[ac)ac cosz- (25) 

^ cos e Mp -I- m 



Assuming that at ac,out, — and the satellite orbital angular momentum is dominated, 
Lzc ^ "^^(ac,out)a„out^ ^ m^/GM^a~^t- From Lzo = Lzc, 



Q \ 2 / N -2 

c,out 2/ 0\/-||<.\2/^\ 2 



~ «M- (l + /or 77- cos^o (26) 



1/3 / _ \ -2 / n \ -2 / /l,f \ -2/3 



Substitution of m = 0.02Mp and 70 = 10° yields a^out ~ 42i?p, which is consistent with the 
numerical results in Fig. [5ll. 



3.4 Classification of tidal evolution 

[Figure [T3] 

Summarizing the above analytical arguments on the tidal evolution, the m/Mp-70 plane is 
partitioned into three regions labeled A, B and C as shown in Fig. [121 In region A, the satellite 
monotonically recedes from the planet until a reaches a^out- With parameters in region Al, 
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a 



oo 

'CjOUt 



is larger than acrit, so that synchronous state = n with 7 ~ i ~ e 



is reahzed. 



With parameters in region A2, synchronous state is achieved before a reaches ~ Ocrit- Then, 
the obhquity and inchnations (7, i and e) do not vary from initial values. In this case, since 
7q ~ 7 ~ -i^ cicout is given by Eq. (j^Tj) without the factor cos^7o. Hence the boundary is 
independent of 70. In region A3, because of small satellite mass, the obliquity evolution is 
dominated by the stellar tidal torque. In regions B and C, the satellite's orbit expands until 
Ccrit- At Ocrit, ^ begins to oscillate by the character change of precessional motion so that the 
satellites begin to migrate inward. The satellite stops at a^out for the parameters in region C, 
while it falls onto the planet in region B. 

Regions B and C are separated by Lzo = 1.755, which is represented by the dashed line in 
Fig. [131 The boundary of Al from B or C is given by acrit = ci'^outy represented by the solid 
line in Fig. [131 If m exceeds the mass determined by the above equation with 70 = 0°, a'^^^^^ 
is smaller than Ocrit for any 70. This separates A2 from B or C. The stellar tidal effects are 
stronger than satellite's one if the timescale (AtAi) for orbital expansion to a'^^^^ is longer than 
the timescale (At as) required for the synchronism by the stellar tide (A^ai and At as are given 
in section 3.5). This separates A3 from Al. 

In Fig. [HI we compare these analytical boundaries with the numerical results in Fig. [91 
Circles and filled squares represent Al and A3 evolution, respectively. In the numerical calcu- 
lation, A3 evolution with 70 < 90° is similar to Al evolution. Triangles, crosses and filled circles 
represents A2, B and C evolution, respectively. The analytically estimated boundaries are in 
good agreement with the numerical results, especially the boundary of Al from B or C and 
that separating B and C. Since the transition of precession from type I to II is not determined 
exactly by Eq. ([7]), the boundary which separates A2 from B or C is not clear enough in the 
numerical calculation. 



3.5 Variation magnitude and timescales of obliquity evolution 

So far we have been concerned with diversity of tidal evolution and intrinsic dynamics that 
regulates the diversity. In this subsection, we discuss ranges of the obliquity changes and their 
timescales that may be important for implications for planetary climate. Typical obliquity 
changes in individual types of evolution are illustrated in Fig. [T3 In general, the variation 
A7 in the obliquity is large for large sin 70. For B and C evolution, not only the precession- 
averaged A7 but also the oscillation amplitude of 7 during precession is large after the satellite's 
migration turns back. 



The evolution timescales of the obliquity in Al and A2 are comparable to migration timescale 
to ac,out- The migration timescale depends on the specific dissipation function Q, which is de- 
fined as the inverse of the frictional energy dissipation per cycle of the tidal oscillation, related 
to the phase shift 5 as l/Q ~ 25 ~ 26t{Q - n) (e.g., Murray and Dermott 1999). With Q, 
Eq. ([TTj) is approximately written as 



[Figure [14] 



[Figure [U] 




(28) 
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Integrating this equation from oq to ac.out is 



13\12nGj k2p^/^ \ m J ]\Rp J V^p. 

(The term of Oq is neghgible.) Substituting Eq. (^71) into Eq. (12^ . the timescale for Al evolution 
is 

Vl0yv0-3y V0.33/ V5-5gcm-3y 

X (1 + /o)^' (^;^) cos'' 70 [year]. (30) 

Timescale for A2 evolution is given by Eq. fl30l) with 70 = 0° for any 70. Figure [16] shows Al 
evolution time calculated by Eq. fl29|) with Q = 12 (current Earth's value) and ^2 = 0.30. Note 
that A^Ai strongly depends on the mass ratio m/Mp and the initial obliquity 70. Also note 
that Q could be larger than the current Earth's value for planets in early stage, so that the 
estimate of evolution timescale includes some uncertainty. 

[Figure [T6] 

The satellite's orbital evolution of type B and C is expansion to a^rit followed by inward 
migration to cic.out or to the planetary surface. The characteristic evolution timescales are given 
by integrating Eq. ( l28l) to Ocrit and multiplying it by a factor 2. Replacing ac by acrit in Eq. (l29l) 
and x2, and substituting Eq. (JTj) into it, 



B,C 



AtRo ~ 1.5x10^ 



loy Vo.3y Vi^0/ viau7 

The A3 evolution timescale is the time necessary for Vt to be decreased to orbital mean 
motion by the stellar torque. Neglecting the first term in r.h.s. of Eq. (ITOl) . 

dVL Sko f MA fR ^ ^ 



dt 2aQ \MpJ \ Op 



' ' nl. (32) 



Integrating Eq. ([32]) from fio to rip = {GM^/aiy/^, 

A*A3 - ^?f^)^^(no-"p) (33) 



3 X 10" {^ \ 



10/ V0.3/ VlM^ / VIAU/ 



(0I3) (5:5^^) (sSi) 
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4 Discussion 



Here we comment on the spin-orbit resonances along the tidal evolution. Other planets' gravi- 
tational perturbations cause precession of the planetary orbit about the invariant plane of the 
planetary system. When the spin axis precession frequency (us) is commensurate with one of 
eigenvalues of the orbital precession frequency (z/q), the spin axis can fluctuate with large ampli- 
tudes. Ward (1973) pointed out that Mars' obliquity has suffered from a large scale oscillation 
because of these spin-orbit resonances. Laskar and Robutel (1993) showed that in the absence 
of massive satellites all of the terrestrial planets could have experienced large-scale obliquity 
variations. Atobe et al. (2004) showed that terrestrial planets in habitable zones in extrasolar 
planetary systems with a gas giant (s) generally tend to undergo the spin-orbit resonance if 
they do not have satellites. They assumed that the planets without satellites underwent nearly 
head-on collisions and had relatively slow spin {D > 20 hours). Due to the relatively large mass 
of the Moon, Earth's spin axis precesses rapidly and avoids the orbital precession frequencies 
of the planets (Ward 1974; Laskar and Robutel 1993). 

As the orbit expands due to the tidal evolution, satellite's contribution to the spin axis 
precession gradually decreases. Initially, z/g ^ z/q ~ lO^'^-lO"^ rad/year, because of small D 
and a. Because a satellite was formed, the impact must have been grazing one, so that initial 
D may be ~ 5 hours or less. Since D and a increase through the tidal evolution, decreases. 
If a satellite is more massive than the lunar mass (~ O.OIM^), the expansion of a is limited 
because of small a^out (Eq- EZj). Then the decrease of stops before it goes down to the level 
of z/q, so that the spin-orbit resonance is avoided. If a sateUite is light enough, the expansion is 
so slow that spin-orbit resonance is not realized during main-sequence lifetime (~ 10^° years) 
of solar-type stars. If a satellite has comparable mass to the lunar mass, Vs eventually becomes 
as small as z/q within 10^° years through tidal evolution, and then obliquity may suffer large 
variation by the spin-orbit resonance. Hence, the spin-orbit resonance may be important for 
planets without satellites (always) and for planets with satellites of about the lunar mass (at 
some point in tidal evolution within main-sequence lifetime). 

So far, we have been concerned with planets at ~ lAU, that is, planets in habitable zones 
around solar-type stars {M^, ~ IMq). Habitable zones around lower-mass stars (e.g., M dwarfs) 
are well inside lAU, because of the low stellar luminosity. As shown in Eq. 0341) . small Op leads 
to rapid removal of angular momentum from the planet-satellite system. Whatever divergent 
tidal evolution the system has undergone, the satellite eventually falls to the planet and the 
planet spin eventually becomes straight and synchronous with its orbital motion by the stellar 
tidal effects within 10^ years. Thus, planets in habitable zones have diversity in planet-satellite 
configurations during main-sequence phase of host stars, if the stars are solar-type stars or more 
massive stars. 

Final accretion stage of terrestrial planets would be multiple collisions of protoplanets (see 
Sec. 1). The probability of a grazing impact to produce a satellite may be higher than head-on 
one. Since the spin produced by the impact is more likely to be tilted, the formed satellite may 
be perished onto the planetary surface through type B tidal evolution on timescales of Myrs 
or trapped near the planetary surface through type C evolution. This means that the satellite 
may orbit near the planetary surface when a next protoplanet approaches the planet. Velocity 
dispersion of protoplanets would be similar to their surface escape velocity, which is also similar 
to Keplerian velocity of the satellite in the proximity of the planetary surface. Hence, the energy 
and angular momentum exchange with the satellite could alter the orbit of the approaching 
protoplanet. This effect might alter scattering/collision cross sections of protoplanets. 
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5 Conclusion 



We have investigated obliquity evolution of terrestrial planets due to tidal interaction with their 
satellites and host stars by numerical integration and analytical arguments, with wide variety 
of initial conditions on the basis of recent N-body simulations of planet accretion. We found 
three domains in the parameters of satellite-to-planet mass ratio (m/Mp) and initial obliquity 
(70) in which the evolution is qualitatively different from one another. 

If a satellite is formed from a debris disk created by an off-center giant impact and the 
planetary spin is dominated by the impact, the planet would have various initial obliquity and 
a satellite with various mass. And planetary spin axis and orbit normal of the satellite are 
almost aligned (e ~ 0). Since N-body simulations show that the impacts are almost isotropic, 
an initial planetary spin axis is more likely to be tilted (sin 70 ~ 1) from the planet orbit 
normal than to be aligned with it (sin 70 ~ 0). In order to study tidal evolution of extrasolar 
terrestrial planet- satellite systems, we numerically integrated evolution with various m/Mp 
and various initial 70, focusing on tilted ones. (We consider the regions at ~ lAU and assume 
eo ~ 0.) Most of previous studies on tidal evolution, however, have focused on a particular 
case of the Earth-Moon systems in which m/Mp = 0.0125 and 70 ~ 10° (e.g., Goldreich 1966; 
Mignard 1978, 1979, 1980; Touma and Wisdom 1994) or the 2-body planar problem 7 = 0° 
(e.g., Counselman 1979). 

As the previous studies show, the satellite orbit first expands with the almost constant 
obliquity (7), the inclination of the satellite orbit normal to the planet one (i), and that to the 
planetary spin axis (e) until orbital radius of the satellite (a) reaches ~ Ocrit ~ 15-Rp (Eq. [7]) at 
which the precession about total angular momentum vector is transformed into that about the 
planet orbital normal. We have found that this character change in precession causes diversity 
of tidal evolution as follows: 

1. At a ~ cicrit, e starts oscillation between ~ and ~ 270. The enhanced e instantaneously 
reduces outer co-rotation radius (cicout)- For large sin 70 cases, a^out strides across a 
without capture at the synchronism and the outward migration of the satellite orbit is 
reversed to inward one. 

2. The evolution with the migration reversal is further divided into two types. The initial 
Lz determines whether the turned-back satellite falls onto the planet (type B evolution) 
or it is captured at ac,out that has jumped inward to ~ 5-10i?p (type C evolution). 

3. With the analytical arguments, the parameter space of m/Mp and 70 is divided into three 
domains (A, B, and C) of qualitatively different tidal evolution (Fig. [T3l) . which is in good 
agreement with the numerical integration (Fig. [H]). Type A evolution is further divided 
into three types: Al) evolution similar to the Earth-Moon system, A2) evolution with 
the synchronism at < Ocrit for m/Mp ^ 0.05, and A3) evolution dominated by stellar tidal 
torque for m/Mp <, 0.005. 

4. In the final state approaching a^out or the planetary surface, e ~ 0°. On the other hand, 
7, i ~ 0° for 7o < 90° while 7, i ~ 180° for 70 > 90° except for A3 evolution. Typical 
7 evolution in each type is shown in Fig. [T3 The variation timescales are evaluated in 
Figs. Uni The timescales for B and C evolution are relatively short (< 10^ years), because 
the satellite orbit turns back at relatively small radius acrit- 

If the formation of a satellite that is comparable to or larger than the Moon's mass and 
planetary spin state are determined by a giant impact, the tidal evolution of the satellite orbit 
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and the spin is most likely to be type B/C that have been discovered by the present paper. (We 
will elsewhere discuss the predicted distributions of m/Mp and 70.) To address dynamics of 
terrestrial planet-satellite systems in extrasolar planetary systems, it is important to consider 
the large diversity of tidal evolution. This diversity may be also important to discuss the 
evolution of surface environments of extrasolar terrestrial planets and their habitability. 
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Appendix A: Equations of precessional motions 

Here we briefly summarize the precessional torques and detailed equations of precessional mo- 
tions derived by Goldreich (1966). The torques acting on the planet (Lp) is sum of those 
from the satellite and the star: Lp = Lps + Lp*. Similarly, torques acting on the satellite is 
Ls = Lsp + Ls*, where subscripts "s", "p" and "*" represent the satellite, the planet and the 
star respectively. Lsp averaged over the orbital periods is given by 

Lsp = - ^^J^ J2Rl{s ■ n)(s X n) = -L(s ■ n)(s x n), (35) 

where a is the semi-major axis of the satellite's orbit, a = {GMp/R^Y^'^, and J2 is the planetary 
oblateness coefficient given by 

- ~2gm; - 2"^' ^^^^ 

{ks is the secular Love number). For the current Earth, ks = 0.95 as to give the known value of 
J2 (Munk and MacDonald 1960). Replacing m and a in Eq. ( l35l) by and Op, the averaged 
precessional torque exerted by the planetary equatorial bulge on the host star (L*p) is 

L,p = _ ^^*^p j^^2(g . k)(s X k) = -Ki{s ■ k)(s X k), (37) 

Op 

where Op is the semi-major axis of the planet's orbit. The averaged torque exerted by the 
satellite on the star (L*s) is 

L,s = -^^^a2(n • k)(n X k) = -K^in ■ k)(n x k). (38) 
4 al 

Since the satellite is regarded as a ring with radius a and mass m in the orbit averaging, Mpi?p J2 
in Eq. (|37|) is replaced by ~ ma^ in Eq. fl38l) . 
Since Lp^ L^p, Lpg Lgp, Lg^ L^g, 

^ = L(s-n)(sxn) + iri(s-k)(sxk), (39) 

^ = -L(s ■n)(sxn)+ 7^2(11 ■k)(nxk). (40) 
dt 

Taking the dot products of Eqs. fl39|) and fHOj) with s and n, respectively, dH/dt = dh/dt = 0. 
These equations yield Eqs. and (jl]). 

Dotting k into Eqs. ([3]) and (jl]) and forming the combination hHd{s ■ n)/dt, 

(41) 
(42) 

V, (43) 



dx 


L 






H 




dy _ 


L 

——zw 
h 


dt 


dz 




'dt 





H' 



where x = s ■ k = cos 7, y = n • k = cos ^ = s • n = cos e and 

w"^ = [(s X n) ■ k]^ = 1 - - - + 2xyz. (44) 
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From Eqs. (HT]) . P2l) and (H3l) . we obtain the new constants 

= Hx + hy, (45) 
X = i^ix' + 7^22/' + iv^'. (46) 

The conservation of arises because the external (stellar) torques on the planet-satellite system 
lie in the planet's orbit, x is a kind of potential energy. Differentiating Eq. (jH]) with Eq. fHTl) 
to g3]) yields 

dw L . . L , . f K2 Ki \ , , , 

Itt " Jj^^y^ -x)- -j^z{xz -y)+ { -^y - —X 1 [xy - z). (47) 

The numerical methods are as follows. Precessional motions are calculated by determining 
initial data for x, y, and z using Eqs. (H5l) . (H6l) and dB]), and w with w = without any loss 
of generality, and simultaneously integrating x (Eq. HI]), z (Eq. |13]), and w (Eq. |47j) with the 
conservation law for Az (Eq. I45l) to eliminate using a 4-th order Runge-Kutta method. Then 
we average the tidal torques given by Eqs. ( l55l) to ( l66i) in Appendix C over one precessional 
period with the resulting (instantaneous) values of x, y, z, and w. Tidal evolution is calculated 
by integrating the tidal equations given by dH/dt (Eq. |48]), dkz/dt (Eq. [51]), da/dt (Eq. 
and dx/dt (Eq. [52]) with the updated averaged torques and x, y, z in the previous step. With 
the updated H, Az, a, and x, we go back to the calculation of precessional motions. 

This method of integration assumes that the three timescales of orbit, precession and tidal 
evolution are well separated. In some cases, the precessional periods can be comparable to or 
longer than the tidal evolution timescale. In the cases, we also integrate directly Eqs. ([1]) and 
([2]), instead of the hierarchical method. We found that the hierarchical method produces the 
results in good agreement with direct method. 



Appendix B. Equations of tidal evolution 

In the planet-satellite-star system, the principal tidal change is brought about by the frictionally 
retarded tide on the planet by the satellite, which results in a loss of mechanical energy from 
the planet-satellite system and angular momentum transfer from the planet's spin to the orbital 
motion of the satellite. Additional smaller change produced by tide due to the host star results 
in dissipation of energy and angular momentum transfer from the planet-satellite system to 
orbital motions about the star. Assuming a phase-locked spin of the satellite on a circular 
orbit, the tide by the planet on the satellite does not lead to any secular change. Thus, we only 
consider planetary tides induced by the star and the satellite. 

In Appendix B and C, all quantities are precession-averaged ones. For simplicity, notations 
for the averaging are omitted. Dotting Eqs. ([5]) by s and n yields 



Dotting Eqs. ([5]) by k yields 



dH ^ dh ^ 
^ = T„.s.- = T,.„ 

~(U ^ H ' 

dy ^ T, ■ k - yXs ■ n 

dt h 



(48) 

(49) 
(50) 
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From these equations, 

Differentiating Eq. (H^ yields 

dx _ 2Ki 

From Eq. (gH]) with H = Jpfi, 



2K, 



xTp ■ k + — ^y(Ts ■ k + yTs ■ n) 
n 



dVL Tp ■ s 



dt Ip 

where Ip is the planet's moment of inertia, given by aMpR"^. Since h = mno?, Eq. 

da Ts ■ n 
— = 2a — ; — . 

dt h 



(51) 

(52) 

(53) 
8]) yields 

(54) 



APPENDIX C: The averaged tidal torques 

In the present paper, we only consider frictional dissipation due to tidal deformation of the 
planet caused by the satellite and the star. The averaged tidal torques on the satellite (Ts) is 
the sum of the torque by the satellite tide (Tgs) and that by the stellar tide (T^^,). If the spin 
axis is not aligned with the orbit normal of the tide raising body, the spin can carry the tidal 
bulge out of the orbital plane. This out-of-plane bulge can produce in-plane torques on a third 
body. Similarly, the torques on the star (T*) is the sum of T^,^. (due to the stellar tide) and 
T^Ks (due to the satellite tide). The torques on the planet are the opposite of the torques on 
the exterior bodies, Tp = -(T, + T,) = -(T„ + T,, + Ts, + Tss). 

Adopting the constant time lag model by Mignard (1981) and Touma and Wisdom (1994), 
in which the distortion of the planet is delayed from the tide raising potential by the time lag 
5t, the torques are 



Tss ■ s 



Tss n 



6t- — ^ 



3 

-fi sin^ e + 3 cos e(fi cos e — n) 



k2Gm^Rl 
St- — - — ^ [3(ficose-n)] 



6t- — ^ 



-i7(cos 7 — cos i cos e) + 3 cos i{Q cos e — n) 



(55) 

(56) 
(57) 



where k2 is the dimensionless tidal Love number (Munk and MacDonald 1960), and 



6t- 



k.GM?Rl 



* p 



6t 



k2GMlRl 



-VL sin^ 7 + 3 cos 7(fi cos 7 — % 



-i7(cos e — cos 7 cos i) + 3 cos i(f2 cos 7 — np) 



6t- 



k^GMlRl 



[3(fi cos 7 — Hp)] 



(58) 
(59) 
(60) 
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Denoting the longitude of the ascending node of the planet's equatorial plane on its orbital 
plane and that of the satellite's orbital plane on the planet's orbital plane as $ and 



and 



koGmM^Rl 



- sin^ i sin^ 7 cos 2($ — ^) — - sin^ i sin^ 7 
8 8 



3 3 
— - cos 7 sin 7 cos i sin i cos($ — ^) + - sin^ 7 



n 







— cos i sin i sin 7 cos($ — ^) 

4 / \ / 



(61) 
(62) 
(63) 



k2GmM^Rl 

fist- p 



- sin^ i sin^ 7 cos 2($ — ^) sin^ i sin^ 7 

8 8 



3 3 

- cos 7 sin 7 cos i sin i cos(<l> — \&) + - sin^ 7 



3 3 

- cos(<l> — \E') sin i sin 7 cos^ i cos 7 cos i sin^ i 

4 4 



(64) 

(65) 
(66) 
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Figure captions 



Figure [T] The tidal evolution of the Earth- Moon system: (a) the obliquity of the Earth to the 
ecliptic 7, (b) the inclination of the Moon's orbit to the ecliptic (c) that to the Earth's 
equator e. The ranges of oscillation during precession are indicated by gray patches. 

Figure [2] Tidal deformation and phase lag are illustrated. The planet (represented by the 
circle in dashed line) is distorted (oval in thick solid line) by the gravitational force of the 
external body (the small circle in thick solid line). The tidal friction causes a phase shift 
(6) in the response of a planet. The case in which the planetary spin is faster than the 
external body's orbital motion is shown. 

Figure [3] The geometry of the orbital plane of the planet, its equatorial plane, and the orbital 
plane of the satellite is shown. Their unit normal vectors are s, k, and n. The obliquity 
7, the inclinations i and e are indicated. 

Figure |4] The trajectories of the tips of s and n projected onto the orbital plane (X-Y) of the 
planet in the case of the Earth-Moon system, (a) a ~ 5R^, (b) a ~ 15Rq, (c) a ~ 60R^ 
(the present configuration). The solid and dashed lines show trajectories of s and n, 
respectively. Filled squares and triangles indicate the positions of s and n. They move 
on the trajectories in labeled number order. 

Figure [5] The numerical calculated tidal evolution of the system with m = 0.02Mp, 70 = 10°: 
(a) the obliquity 7, (b) the inclination of the satellite orbit to the planetary orbit i, (c) 
that to the planet's equator e, and (d) the semi-major axis of the satellite a (dashed line), 
the inner co-rotation radius a^in (gray solid line) and the outer one ac,out (black solid line) 
calculated by Lzo = Lzc (Eqs. [2^ and 1251) . In the panels (a), (b), and (c), the ranges of 
oscillation during precession are indicated by gray patches. 

Figure M Same as Fig. El but m = O.OlMp, 70 = 80°. 

Figure [7] Schematic illustration of the precession motion of s and n. When a < acrit, s and 
n precess about their total vector (type I precession). During the precession, the mutual 
inclination (e) between s and n is almost constant. When a > Ocrit, s and n precess about 
k, respectively (type III precession). During the precession, e oscillates from 0° to ~ 27. 

Figure M Same as Fig. but m = 0.04Mp, 70 = 40°. 



Figure [9] The dependence of tidal evolution on initial conditions (initial obliquity 70 and 
satellite-planet mass ratio m/Mp) found by the numerical integration. Type Al evolution 
with 7 0° (for 7o < 90°) or 7 ^ 180° (for 70 > 90°) is plotted by open circles. A3 
evolution which shows 7^0° with 70 > 90° is plotted by filled squares. Crosses and 
filled circles represent B and C evolution, respectively. Triangles represent A2 evolution. 

Figure 1101 Contours of angular momentum L (solid lines) and energy E (dashed lines) in a 
planar system. Planetary spin and satellite's orbital periods are equal when Cl = n^, 
represented by the dotted line. Critical points are shown by filled triangles, at which all 
lines of constant-L, constant-i?, and synchronism are mutually tangent. Path of tidal 
evolution in the {Q, n)-plane are indicated by arrows. 

Figure 1111 The time evolution of the normalized total angular momentum (L) and its Z- 
component (Lz), and the total mechanical energy of the planet-satellite system (E), in the 
cases of (a)m = 0.02Mp and 70 = 10°, (b)m = O.OlMp and 70 = 80° and (c)m = 0.04Mp 
and 70 = 40°. 

Figure 1121 Path of the tidal evolution of non-planar systems in the case of (a) m = 0.02Mp 
and 70 = 10°, (b) m = O.OlMp and 70 = 80°, (c) m = 0.04Mp and 70 = 40°. Thick solid 
lines are the numerical results. Solid and dashed lines are contours of L lines assuming 
e = and E. Spin-orbit synchronism {Cl cos e = n^) is represented by the dotted lines. 

Figure 1131 Qualitatively different tidal evolution predicted by the theoretical arguments, as a 
function of (m/Mp,7o). B and C are separated by \Lz\ — L* (thick dashed lines). Al 
and B are separated by a^it — c^^out (thick solid line). A2 and B/C are separated by the 
vertical line which has the mass determined by a^nt — (^Tout with 70 = 0°. The boundary 
between Al and A3 is given by AtAs = AtAi (thin solid lines). 

Figure 1141 Comparison of the analytically estimated boundaries with the numerical results in 
Fig. [9l For the analytical estimate, see text. 

Figure 1151 Schematic illustration of typical obliquity changes in individual evolution types. 

Figure [16] The timescales of the tidal evolution, with Q = 12 and ^2 = 0.30, using Eqs. fl30|) . 
The contours of the timescales are represented by dotted lines. From right to left, contours 
represent 10^, 10^, 10^^ yrs. 



27 




28 



Figure 2: Atobe and Ida 
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Figure 3: Atobe and Ida 
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Figure 4: Atobe and Ida 
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Figure 6: Atobe and Ida 
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Figure 8: Atobe and Ida 
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Figure 9: Atobe and Ida 
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Figure 11: Atobe and Ida 
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Figure 13: Atobe and Ida 
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Figure 14: Atobe and Ida 
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Figure 15: Atobe and Ida 



42 




43 



